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Abstract. We analyze the open boundary partially asymmetric exclusion 
process with smoothly varying internal hopping rates in the infinite-size, mean 
field limit. The mean field equations for particle densities are written in terms 
of Ricatti equations with the steady-state current J as a parameter. These 
equations are solved both analytically and numerically. Upon imposing the 
boundary conditions set by the injection and extraction rates, the currents J 
are found self-consistently. We find a number of cases where analytic solutions 
can be found exactly or approximated. Results for J from asymptotic analyses for 
slowly varying hopping rates agree extremely well with those from extensive Monte 
Carlo simulations, suggesting that mean field currents asymptotically approach 
the exact currents in the hydrodynamic limit, as the hopping rates vary slowly 
over the lattice. If the forward hopping rate is greater than or less than the 
backward hopping rate throughout the entire chain, the three standard steady- 
state phases are preserved. Our analysis reveals the sensitivity of the current to 
the relative phase between the forward and backward hopping rate functions. 
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1. Introduction 

Asymmetric exclusion processes (ASEP) have been used as model nonequilibrium 
statistical mechanical systems to represent many physical processes such as traffic 
flow mGlE], ion transport across channels 01 Ej, mRNA translation [El El IE], and 
vesicle translocation along microtubules For uniform hopping rates, the steady- 
state currents, particle densities, and correlations of the one-dimensional totally 
(TASEP) and partially asymmetric exclusion process (PASEP) have been studied 
extensively using recursion methods, exact matrix product techniques, and mean field 
approximations |H3 EU E2 El 03] The PASEP is described in Figure □ and is 
comprised of one-dimensional lattice of sites, each of which can only be empty of 
singly occupied. The rules governing the dynamics of this system are as follows: a 
particle at site n hops to site n + 1 with probability p n dt in the time infinitesimal dt, 
only if site n + 1 is empty. Similarly, it can hop backward to site n — 1 (if site n — 1 
is empty) with probability q n dt. At the left and right boundaries (sites n = 1 and 
n = N, respectively) the injection probabilities are adt and 8dt, respectively, provided 
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these sites are unoccupied. Extraction of particles from sites n = 1 and n — N occur 
at rates 7 and /?, respectively. Particles do not hop if others are blocking their target 
sites. We will only consider the averaged steady-state configurations of this system. 

Exact steady-state currents for the case of constant p n — p and q n = q have 
been found In the case where p 7^ q, the exact solution in the infinite 

lattice limit (N — > 00) exhibits three phases described by maximal current, high 
particle density, and low particle density. Within each of these phases, the steady- 
state particle current is described by explicit analytical expressions ^2]. Additional 
subphases corresponding to different density profiles arise within the high and low 
density current regimes |15l I16| . When the forward and backward hopping rates (out 
of and into each site) are equal, the chain is purely diffusive and is driven only by a 
difference between the injection/extraction rates at the boundaries. In this case, only 
a single, smooth (with respect to the injection/extraction rates) current phase exists. 

In many systems modelled by the PASEP, the internal hopping rates are spatially 
varying. For example, variations in the hopping rates may arise in pores that have 
internal molecular structure, microtubules tracks (on which molecular motors move) 
that are comprised of periodic subunits, or from variations in mRNA or DNA sequence. 
Variations in the forward hopping rate for fixed lattice defects in a TASEP have been 
treated approximately in the limit of few, isolated defects |H1 117) , and in the periodic 
case where the forward hopping rate takes on two values [T5| . 

In this paper, we consider spatially varying internal forward and backward particle 
hopping rates in a PASEP. We find solutions for the current and density when the 
forward and backward hopping rates are given by functions p n and q n that vary slowly 
with the lattice position n. In the thermodynamic, mean field limit, the equation of 
motion for the mean occupancy at each site can be described in terms of a nonlinear 
continuum equation involving the coarse grained mean occupations cr(x), and the 
continuum hopping rate functions p(x) and q(x). In the next section we derive the 
steady-state continuum equations by expanding the occupancy evolution equations 
in powers of e = l/N, where N — > 00 is the total number of lattice sites in the 
chain. We the consider four general classes of the hopping functions p(x) and q(x). In 
Section 3, we treat the "pure diffusion" limit where p n = q n +i, or, in the continuum 
limit, p(x) = q(x + e). In this limit, q(x) ~ p(x) — ep'(x), the mean-field equations 
become linear, and exact simple results are recovered. In Section 4, we consider the 
"shifted diffusion" limit, where the forward and backward hopping rates at each site 
are identical in the sense that p„ = q n , or, \p(x) — q(x)\ -C o(e). We find exact implicit 
solutions for special forms of p(x),. These results are markedly different from those 
found for pure diffusion, although the structure of p(x) and q(x) are nearly identical 
for the two cases. The case where \p(x) — q(x)\ > 0(e), and p(x) — q(x) does not 
change sign is considered in Section 5. Asymptotic analysis indicates that the standard 
three phase structure found for constant p, q is preserved qualitatively. In 

Section 6, we show that internal density boundary layers arise if p(x) — q(x) crosses 
zero for < x < 1. This case eluded analytic treatment so only numerical and 
simulation results were obtained. In all cases, we compare our results with numerics 
and continuous time Monte-Carlo simulations. In the Summary and Conclusions, we 
discuss the limits in which one would expect mean-field approaches to yield exact 
steady-state currents. 
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Figure 1. The jV+l-site open boundary, partially asymmetric exclusion process. 
The continuum limit is a taken by setting each lattice site to size e = 1/N, thereby 
normalizing normalizing the total length. 



2. Continuum Mean Field Limits 

Consider a one-dimensional lattice (Fig. 1) containing N + l sites each of length e. For 
the interior sites, the continuum limit of this lattice will be defined by a sampling of all 
relevant quantities (e.g., density) at the centers of each lattice site. Density profiles in 
the presence of sources and sinks exhibit rich shock behavior as studied by Parmeggiani 
et al. [T^l, and by Evans et al. [201 ■ Here, we will neglect adsorption/desorption at 
the interior sites; however, we allow the internal hopping rates to vary slowly along 
the chain. 

The equation for the discrete occupation variable a n € (0, 1) in the chain interior 

is 

= (ft-1 - ft) + (ft+1 ~ Jnl l<n<JV-l, (1) 

where 

J+ = p„<7 n (l - & n +i) and J~ = q n fr n (l - & n -i) (2) 

are the currents from site n to site n+l and from site n to site n — 1, respectively. 

The mean field assumption implies that the ensemble averaged occupancies are 
uncorrelated, (a n a m ) w a n a m , where a n = (<7„). Upon taking an ensemble average 
of Eq. ^ an d applying the mean field approximation, the evolution equation for the 
mean occupancy in the chain interior (l<n<iV — 1) becomes 



^ — ift-l ft) + (J-n+l Jn ) 



(3) 



where J+ = (J+) = p n cr n (l - <r n+1 ) and J n = ( J n ) = q n a n (l - cr„_i). Upon 
extrapolating the continuous function according to a(x = ne) = a n , and Taylor 
expanding Eq. [3] in powers of e, wc find the continuum mean field equation: 



da(x) 



= e[(q- p)a(l a)]' + £ - [[(p + q)a]' (1 - a) + (p + q)a*>] ' . (4) 
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Assuming the steady-state limit, and integrating the conservation Eq. we obtain 

EU 

(p - q)a(l - a) - | ([(p + q)a]' (1 - <r) + (p + q)aa') = J, (5) 

where the integration constant J is the steady-state current. Equation can be 
rewritten in the Riccati form 



ea'(x) = -JP{x) + Q(x)a{x)(l - a(x)), 



where 



P(x) 



Q(x) 



(p + q) 
2(p-g) 



and 



(p + 9)' 



+ (?) (p + 9) 



(6) 



(7) 



The boundary densities are found by measuring the steady-state current into and out 
of the first and last sites: J = a(l — Co) — 7°o and J = (3a n — 5(1 — <tn), from which 
we find 



a — J 



and <7/v 



J 



(8) 



a + 7 " + <5 

Equations and form the basis of our steady-state analysis. Integrating Eq. 
El from x = to x = 1, and imposing the boundary conditions (Eqs. implicitly 
determines J. Once the steady-state current J is fixed, the mean field density profiles 
are determined. For certain p(x),q(x), one may be able to solve Eq. analytically, 
and use this result along with the boundary conditions [S] to find J in closed form. 



3. Pure diffusion: p(x) = q(x + e), < x < 1 

Consider the special case p n = q n +i, 1 < n < N — 1 where the hopping rates between 
two sites are equal. In this case, there is no driving force on the particles and a net 
current arises only from differences in injection and extraction rates at the two ends. 
The quadratic terms in Eq. ^ cancel, the equation for a n becomes linear, and the 
mean field approximation is exact. In the continuum approximation p(x) = q(x + e), 



P(x) 



1 



£ p 



1 + 0- + -r ^ 

2 p 4 



p{x) 

and to lowest order in e, 
a'(x) ~ 



P' 



P_ 

P 



and Q(x) = 0(£ 3 ), 



ep(x) 



Integration of Eq. ^| yields 
<r(l)-<r(0) 



-dx' 



J 



1/P>- 



P(x') 

Upon applying Eqs. |H1 (with q\ = po and = PN-l), and solving for J, 

(a/3-7,5) 



J 



<l/p)(a + 7)(/3 + <J) + {a + p + ^ + S)' 



(9) 



(10) 



(11) 



(12) 
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This same "homogcnization" result, with the Ricmann equivalent 

JV-l 

(l/p)= J2Pn\ (13) 

n=a 

is also easily obtained by recursively solving the exact discrete equation p„_ i(er„ — 
cr n _i) = — J, (1 < n < N). The corresponding density profile is obtained through 

n-1 

<T n =a Q + jJ2 P j 1 . (14) 

For constant p n = q n = p, Eq. I12l rcduccs to trivial result for the steady-state current 
of an N + 1 site boundary-driven chain 

Ar(a + 7 )(/3 + (5)+p(a + /3 + 7 + 5)' 1 j 

4. Shifted diffusion: p{x) ~ </(a;), < x < 1 

If p n ~ g„, such that \p(x) — <C o(e), particles at each site hop equally to the 
right or the left, with possibly different rates from site to site. This case corresponds 
to particles that cannot distinguish forward from backward motion but, as detailed 
balance is violated, is not equivalent to pure diffusion. We will see that this slight 
change in the hopping rate structure from the p n = q n +i case results in a very different 
steady-state current. 

To lowest order, when p(x) ~ q{x), 

P(z)~J- and Q(x)^~e^l. (16) 

pyx) p{x) 

Since Q(x) is of order P(x), the cr(l — a) term in Eq. cannot be neglected and, 
unlike the purely diffusive case, the problem is nonlinear. Therefore, we would not 
expect the mean field particle densities or currents to be necessarily exact. 

In this case, there are various variable transforms that render the Riccati equation 
analytically tractible. The simplest case is where the Riccati equation is separable. 
This occurs when (p + q) 1 —constant, which implies p(x) — q(x) = ax + 6. Equation 
EJcan then be integrated from x = and a = cr(0) = ero to give 

da;' 1 f 7 ^ der 

T. ( 17 ) 



where a± = 1/2 ± 1/2^1 + 4J/(ea). Integrating (|17J) . we find the density profile a(x) 
from 

/ a(x)-a+ \ / (t -(t_ \ 
a(x) - ct_ J \a - a + J ' 

An implicit formula for J (to order e) is found by imposing the boundary condition 
at x = 1 (er(l) = <j N w 5/(fl + 5)): 



+ b\^ 1+AJ/{ea) _ f 5-((3 + 5)a + \ f g-{a + i) a J 



(19) 



The solution to Eq. ^5]is found numerically and plotted in Fig. for representative 
parameters. 

Next, consider another analytic solution found by using the definition 
e y'(x) 

which transforms Eq. to 

Q(x)y"(x) - [Q'(x) + e~ 1 Q 2 (x)] y'(x) + e~ 2 JP(x)Q 2 (x)y(x) = 0. (21) 
Provided 

Q'(x)+e- 1 Q 2 {x) =0, (22) 
and Q(x) ^ 0, we find 

y"(x)+s- 2 JP(x)Q(x)y = 0. (23) 
The condition Eq. [521 is solved by 

W X + b (p + q) (p + q) 



which constrains q(x) to p{x) through 

e -2x/e 



q(x) 



X 

where 



(x' + b)g(x')e 2x '/ e dx' + constant 



(25) 



g(x)=(*-^ n ^p(x)-p>(x). (26) 

Given pairs of p(x), q(x) that satisfy Eq. [2]or[2Sl one can find analytic solutions to Eq. 
1231 reconstruct a(x) via Eq. [201 an d impose the boundary conditions on cr(0) and <r(l) 
(Eqs. [HJl to find an implicit equation for J. Note that the constraint Eq. 1251 allows for 
analytic solutions of Eq. Elfor hopping rate functions more general than p{x) — q(x). 
Restricting ourselves to p(x) = q(x), the only solution for p(x) + q(x) = 2p(x) that 
satisfies Eq. Oil is 

p(x) =a/(x + b). (27) 

Equation 1231 then becomes 

y"(x) + k 2 y{x) = Q 1 k 2 = - (28) 

ea 

admitting a solution of the form 

y(x) oc e lkx +ce~ ikx . (29) 
Upon setting a(0) = Q- 1 (0)(2/'(0)/y(0)) = a = (a- J)/(a + y), 
kb + iao 



kb — iao 

Substituting Eq. OH into Eqs. [Hand [211 we find 



(30) 



, , , en cos kx — kb sin kx 

a(x) = (x + b)k-^-— — . (31) 

do sin kx + kb cos kx 




Figure 2. Densities and currents from analytic solutions of the Riccati equation 
and from simulations, (a) The density profile resulting from the linear hopping 
rate function p(x) = q(x) = x + 1. Results from the analytic solution Eq. 1181 f solid 
curve) and from Monte-Carlo simulations (circles) are shown. Also shown is the 
exact density for the purely diffusive case p(x) = q(x + e) = x + 1 (dashed curve), 
(b) Currents derived from both the analytic solution Eq. 1191 (solid curve) and 
Monte-Carlo simulations (circles). Also shown for contrast is the exact current 
in the purely diffusive case (dashed curve), (c) The density profiles associated 
with the hopping rate function p(x) = q(x) = l/(x + 1). The solid, circled, 
and dashed curves correspond to analytic, Monte-Carlo, and exact diffusion (for 
p(x) = q(x + e) = l/(x + l)) solutions, (d) Currents derived from both simulation 
and the analytic solution to Eq. 1321 The arrows in (b) and (d) mark the value 
a = 1.5 used in plotting the density profiles shown in (a) and (c). 



Finally, imposing the boundary condition at x = 1 implicitly determines J: 

(1) k(b + 1) a ° C ° S ^ _ ^ S * n ^ 5 + J /^2) 

do sin k + kb cos k (3 + S 

Since J ~ ea, Co ~ otj (a + 7) and ct/v ~ S/(J3+S) can be used to numerically solve Eq. 
I32lfbr currents and densities. Expanding J»0 also shows that J oc ea{a(3{b+l)— try 8). 

In Figs. |2J we plot (a) the densities and (b) the currents for the hopping rate 
profile p(x) = q(x) = ax + b as a function of driving a. The results of extensive 
Monte-Carlo simulations using the BKL continuous time algorithm [22 , for a lattice 
of size N = 1000, are also shown. Both analytic (in the e — ► limit) and simulation 
results agree to a high degree of accuracy. In Figs. |2fc) and (d), we plot the densities 
and currents corresponding to the inverse hopping rate profile p(x) = q(x) = a/(x + b). 
Here, the current also agrees well with the Monte-Carlo simulations. However, there 
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is a small discrepancy between the densities from mean field theory and those from 
Monte-Carlo simulations. This discrepancy is not unexpected since correlations are 
neglected in mean field theory. Also shown for comparison are the densities and 
currents for the purely diffusive chain where p(x) = q(x + e). These densities (dashed 
curves) are very close to those corresponding to p(x) = q(x), however, the diffusive 
currents are significantly different. By shifting the backward hopping rate function by 
e, the density profile changes only slightly. However, since the steady-state currents 
scale as e, small changes in boundary densities can lead to large relative differences in 
the steady-state currents. 



5. Completely driven chain: \p{x) — q(x)\ ^ O(e), < x < 1 

In this Section, we consider significantly different forward and backward hopping rates, 
and for simplicity, first assume p(x) > q(x) for < x < I. In this case, neither P{x) 
nor Q(x) is small, but Eq. [5] can be treated using singular perturbation theory and 
the appropriate implementation of density boundary layers. Suppose a boundary layer 
arises near x ~ s. Rescaling x — ey, we find 

= ~JP( £y ) + Q(ey)a(y)(l a{y)). (33) 

ay 

Within the boundary layer, y ~ 0(1 ), P(ey) w P(0), and Q(ey) w Q(0). Equation 1551 
can be integrated to find the left inner solution 

m a + (0)(a o MO)) - q_(0)K - cr + (0))e^- W- + W)QWv 

e [V> ao- ( 7_(0)-(ao-a + (0))e(^( )- ff +( ))0( ^ ' [ ' 

where a± (0) is also the outer solution to Eq. 




evaluated as x — > 0. Of the two possible outer solutions, only cr + (0) can match the 
inner solution <7" l (y —* oo). The uniform solution with a density boundary layer at 
x ~ e is thus 

a e (x) = ap(x/s) + <7+(x) -a + (x = 0). (36) 

This solution automatically satisfies the boundary condition at x = 0: erg(0) = era- 
The current is determined by satisfying the boundary condition at x = 1; and since 
af{y = 1) ~ a+(0), 



P + 5 

The only possible solution to Eq. 1371 is 




(37) 



J= \{ p - 5 -W) [p+5)2 ) + U^- 5 -W) {p+8? 
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In addition to this result, two other solutions to Eq. 1331 exist. One with a 
boundary layer at x ~ 1, and another with boundary layers at both x ~ e and x ~ 1. 
If a boundary layer exists only at cc ~ 1, only <7-(x) can match the inner solution near 
x = 1 and the uniform the solution analogous to Eq. 1361 is 

a r {x) = al n (x/s)+a^(x)-a-(x = 1), (39) 

where 

in, M _ Mi)(^ - Mi)) - Mi)(^v - (r+ (i)) e (^(i)-^(i))q(i)(i-^)/ e 

In this case, the self-consistent current is found from oy(0) = oq: 



J= K a_7 "SI (a+7)2 ) + W(^ 7_ ^i (a+7)2 ) +4aT (41) 

When both boundary layers exist, the outer solutions must match at at least one 
intermediate interior position a+(x*) = a—(x*), <C x* <C 1. The corresponding 
uniform solution is 

( oi(x) < x < x* 

tr.(») = | (42) 
^ ay (a;) x* < x < 1 

with corresponding maximal current 

J - (43) 

Jmax ~ 4P(x*y [ ' 

The solutions Eqs. |2H1^2 and 1431 are valid only in the parameter regimes where 
J < Q(x*)/4P(x*) and < a < 1, and represent a generalization of the well- 
established three phase current structure arising in the constant p, q PASEP ^| • 
This three phase structure is preserved only if p(x) — q(x) does not change sign on 
x £ [0, 1]. Note that if p, q are constant and p — q > 0, the inner solutions are exact 
on x £ [0, 1] and we recover the known results for the PASEP ^3^]. The mean field 
densities and currents for p(x) —2 + 2(x — 1/2) 2 and q = 1/2 are plotted in Figurc|21 
along with results from continuous-time Monte-Carlo simulations. The agreement is 
extremely good between the asymptotic mean field currents and simulation currents, 
suggesting that the basic physics of the three phase structure is preserved and that 
mean field theory provides exact, steady-state currents. The densities are also in good 
agreement, except in barely discernible region within the boundary layers where mean 
field and simulation derived densities differ. Note that there is also a slight discrepancy 
between mean field and simulation currents in the maximal current regime (Fig. 3c). 
The underestimation of the current by the mean field analysis results from the finite 
size of the rate-limiting region. Although s is small, and p(x) is reasonably slowly 
varying, the rate limiting region at x w 1/2 is small enough for actual current (from 
MC simulations) to be noticeably greater than the asymptotic mean field result. 
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Figure 3. Densities and currents from asymptotic solutions of the Riccati 
equation and from simulations. Parameters used are e = 1/1000, j3 = 7 = 1, 
8 = 0.5, and p(x) = 2 + 2(x - 1/2) 2 > q(x) = 1/2. In all plots, the solid 
black curves correspond to asymptotic solutions of the mean field equations, while 
the blue circles correspond to results provided by Monte-Carlo simulations, (a) 
Asymptotic and simulation densities for a = 0.5. These parameters used render 
the system in a low density, entry rate-limited regime, (b) For a = 1.5, the system 
is in the maximal current phase (J m ax = 3/8), where boundary layers arise at 
both x fa and x ss 1, and x* = 1/2. (c). The steady-state current as a function 
of a. Given the other parameters used, the system transitions from a low density 
to a maximal current phase at a = 5/6. The parameters a = 0.5,1.5 used in 
plotting the densities in (a) and (b) are marked with arrows. 



6. Opposing drifts: \p(x) — q(x)\ 0(e) except at countable points x$ 

Finally, consider the important class of hopping rates where \p(x) — q(x)\ 3> 0(e), 
except at certain points xq where p(x) — q(x) crosses zero, like Q(x) in the e —> limit. 
Examples of p(x), q(x) with these properties are p(x) = a + bx, q(x) = a + 6(1 — x) 
(where x = 1/2), and periodic p(x),q(x) such that Q(x) oscillates above and below 
zero. Periodic hopping rates may arise during transport through pores with atomic 
periodicity. For example, periodic arrangements of atoms or molecules within the 
pore would impart a periodic potential on translocation of particles of the form 
p(x) cx exp[(V(a;) - V(x + e))/k B T] and q(x) cx exp[(V(x) - V(x - e))/k B T], which 
are periodic if V(x), the interaction potential as a function of the coordinate x along 
the axis of the chain is itself periodic. 

Instead of specifying a detailed, molecular model for the hopping rates, we assume 
p(x),q(x) to be functions that qualitatively capture the physics arising from periodic 
hopping rates. The qualitative dependence of the steady-state currents and densities 
should not depend upon the exact, quantitative forms chosen for the periodic hopping 
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Figure 4. Densities and currents from solutions of the Riccati equation and from 
simulations for the periodic hopping rates modeled by Eq. 1441 The parameters 
used are: a = 1, 6 = 0.5, a = 1.5, /9 = 7 = 1, 5 = 0.5. (a) p(x) (dashed), 
q(x) (dotted), and the numerically solved density profile for k = I0,<f> = 0. 
(b) Numerically computed density profiles for k = 11.3, <j> = (solid) and 
k = 11.7, <j> = (dashed), (c) Numerically solved (solid) and simulation-derived 
(circles) steady-state currents as a function of k (tj> = 0). The arrows indicate 
k = 11.3, 11.7 used for generating the profiles in (b). Currents associated with 
integer values of k are denoted by filled circles, (d) The current as a function of 
the phase difference <f> between p(x) and q(x) for k = 12. 



rates. Therefore, for simplicity, we assume: 
p(x) = a + 6sin 2 (7rfcx) 

q(x) — a + b cos 2 (nkx + Tr(j>), 



(44) 



with k > 1 and a > b. This functional form for the hopping rates captures the 
periodicity of the pore potential and allows for a phase difference 4> between the forward 
and backward hopping rates. When = 0, p(x)+q(x) = 2a+b, P(x) = 2/ (2a + b), and 
the function Q(x) = — 26cos(27rfcx)/(2a-|-&) crosses zero at points x = (2n + l)/(4fc). 
Despite this simplification, there is no analytic solution to Eq. and we were unable 
to find approximations. Asymptotic analysis of the Riccati is also difficult due to the 
existence of multiple, interior boundary layers, and the fact that we must determine 
the boundary densities to 0(e) in order to extract the current J. Moreover, numerical 
solutions to Eq. arc difficult to obtain for extremely small £ since numerical 
errors build up as one integrates Eq. from x = to x = 1. Nonetheless, we 
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compare currents and densities derived from numerics and continuous-time Monte- 
Carlo simulations. 

Figures 4 show numerically computed (e = 0.001) and simulated (N = 1001) 
densities and currents for periodic hopping rates Eq. 1441 In Fig 4(a) are the functions 
p(x), q(x), and the density profile for k = 10 and = 0. Due to the oscillatory nature 
of p(x) and q(x), the densities are locally compressed and rarefied, rapidly jumping 
between a(x) = and u(x) = 1. The numerically computed densities for noninteger 
k are also shown in Fig. 4(b). Small changes in k can cause large variations in the 
density near the x — 1 boundary, causing dramatic changes in the current, as shown in 
Fig. 4(c). Fig 4(d) shows the sensitivity of the current to variations in the phase cf>. For 
clarity have shown only numerically computed density profiles: in our plots, densities 
found from simulations are nearly indiscernable from those found numerically. As in 
all other cases of slowly varying hopping rates, mean-field theory appears to yield 
exact steady-state currents. The agreement between numerical and simulated data is 
extremely good in Fig. 4(c), but the discrepancy increases as the number of hopping 
rate oscillations k increases. 

7. Summary and Conclusions 

We have formulated the mean field approximation of a partially asymmetric exclusion 
process in the hydrodynamic limit in terms of the solution to the Riccati equation. 
This nonlinear equation can be solved in special cases and asymptotically analyzed 
in others. We compare numerical and analytical results with results from extensive, 
(10 9 — 10 1 steps) continuous time Monte-Carlo simulations and find extremely good 
agreement for the steady-state particle currents. The numerical simulations fall well 
within the within the simulation error, typically < 1% and barely discernible in our 
plots. This agreement holds for all hopping rate profiles considered, provided they 
do not vary rapidly along the lattice. Moreover, although we have not proven that 
soutions to the Ricatti equations (Eqs. EJ yield exact currents, comparison of the 
numerical solutions for the current with those obtained from extensive continuous-time 
MC simulations shows a decreasing discrepancy as e/£ — » 0, provided sufficient long 
simulations are performed. Therefore, we conjecture that mean field approximations 
provide asymptotically exact steady-state currents as long as the hopping rate structure 
is smoothly varying in the thermodynamic (N — > oo) limit. 

The simulated densities, as expected, arc quantitatively different from those 
obtained from the numeric or analytic solution of the Ricatti equation. Moreover, 
we find that the three-phase current structure of the PASEP is preserved when 
p(x) > q(x). For cases where \p(x) — q(x)\ < 0(e) (pure diffusion and shifted diffusion) , 
J ~ e, but is sensitive to even slight shifts between the functions p(x) and q(x). The 
cases p(x) = q(x) and p(x) = q(x + e) correspond to physically realizcablc systems, 
yet yield very different results. When p(x) and q(x) vary periodically, as might be 
expected along a molecular channel constructed from a periodic array of atoms or 
molecules, the currents derived from solving Eq. El also appear to be exact, provided 
there are a large number of sites in each period. 
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